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MAXIMUM-LIKELIHOOD  ESTIMATION 


OF  TIME-VARYING  DELAY 


1.  INTRODUCTION 


In  their  classic  paper,  "The  Generalized  Correlation  Method  for 
Estimation  of  Time  Delay,"1  Knapp  and  Carter  presented  the  solution  to  the 
problem  of  maximum-likelihood  (ML)  estimation  of  constant  delay,  dQ, 
between  signals  received  at  two  spatially  separated  sensors  in  the  presence 
of  uncorrelated  noise.  The  received  waveforms  were  modeled  mathematically  as 


ra(t)  =  s(t)  +•  na(t),  -«<  t  <  +«  , 
rb(t)  =  as(t  -  d0)  +•  nb(t),  -  •  <  t  <  +  «  , 


(1-la) 

(1-lb) 


where  a”  was  a  relative  attenuation  constant  and  s(t),  n  (t),  and  n.(t) 

a  D 

were  uncorrelated ,  stationary  Gaussian  random  processes.  Knapp  and 

Carter1  showed  that  the  ML  estimator  of  dQ  can  be  realized  by  a  pair  of 

prefilters  followed  by  a  crosscorrelator.  Their  solution  is  identical  to 

2  3 

that  proposed  by  Hannon  and  Thomson  ’  whose  motivation  for  estimating 
delay  was  to  improve  estimates  of  the  spectra,  cross  spectra,  and  coherence 
of  stationary  time  series.  Here  we  present  a  major  generalization  of  these 
preceding  analyses  that  includes  (1)  arbitrary  time-varying  delay,  d(t); 
(2)  nonstationary  random  signal  process;  and  (3)  arbitrary  observation 
interval . 


Previous  attempts  to  extend  the  theoretical  solution  described  in 
references  1,  2,  and  3  have  been  relatively  limited  in  scope.  Assuming  a 
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stationary  signal  process  and  an  infinite  observation  interval,  Knapp  and 

4 

Carter  obtained  an  approximate  ML  estimator  of  both  differential  delay 
and  Doppler  for  a  source  whose  relative  velocity  is  much  less  than  the 
signal  propagation  velocity.  Here  the  estimator  structure  included  a 
time-compander  following  one  of  the  two  prefilters  before  crosscorrelation 


to  compensate 

for  the 

Doppler  time  scaling  of 

the 

waveform. 

Wax5 

generalized  the 

analysis 

to 

include  differential 

phase 

Several 

other 

authors6'9  have 

described 

the 

degradation 

in  compensated 

and  uncompensated 

crosscorrelator 

outputs 

due 

to  motion 

of  the 

source 

Beyond 

these 

relatively  limited  theoretical  studies,  there  has  been  both  a  need  for  and  a 
continuing  effort  to  develop  practical  algorithms  that  estimate  time-varying 
delay  more  generally.11^  16  The  new  and  general  theory  presented  here  can 
provide  guidance  to  that  effort,  as  well  as  fresh  insights  into  previous 
theoretical  results. 


We  model  the  problem  of  time-varyi ng-delay  estimation  as  follows: 
A  vector  of  real  waveforms, 


r(t)  = 


n't) 

s(t) 

W]  ( t ) 

= 

■t- 

r2(t) 

—  — 

as(t  -  d(t) ) 

_ 

w2(t) 

_  X 

(1-2) 


is  observed  on  the  interval  [ ,  T^J,  where  T.  and  denote  initial 
and  final  observation  time,  respectively.  For  convenience,  we  define  r(t) 
as  zero  for  t  outside  this  interval.  The  signal  s(t)  is  a  sample  function 
of  a  zero-mean  Gaussian  random  process  having  covariance  function 


Rs(t1.t2)  =  E  {  s  ( t  •)  )  s  ( 1 2 )  } 


(1  -3) 
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The  delayed  and  attenuated  signal  as(t  -  d ( t ) )  is  related  to  s(t) 
through  a  nonrandom  but  unknown  invertible  linear  operator, 

Ld(t),a{s(t)>=  as(t  -  d(t) )  .  (1-4) 

The  noise  waveforms  w^(t)  and  w^(t)  are  sample  functions  of  white 
Gaussian  random  processes  having  covariance  functions 


R  (t,.t,)  =  R 


w 


w2(trt2)  -  ^(t,  -  t2) 


(1-5) 


The  signal  process  and  noise  processes  are  mutually  independent .  The 
attenuation  factor  a  and  delay  function  d(t)  in  (1-2)  and  (1-4)  are 
nonrandom  but  unknown.  Since  d(t)  represents  delay,  we  will  assume  here 
that  d(t)  >  0.  This  restriction  can  be  removed  with  a  somewhat  more  lengthy 
analysis.  The  attenuation  constant  ”a  can  be  any  nonzero  real  number.  The 
problem  is  to  estimate  d(t)  and  a*. 


The  model  in  (1-2)  through  (1-5)  assumes  white  noise  processes  and  a 

nond i spers i ve  (frequency  independent)  propagation  medium  for  simplicity. 

One  can  extend  the  developments  of  the  following  sections  to  include 

nonwhite  noise  by  applying  noise  whitening  techniques  similar  to  those 

described  in  reference  17,  p.  290.  One  can  also  include  dispersion,  as  well 

as  time  varying  delay,  by  replacing  the  operation  (1-4)  by  a  more  general 

1  8 

invertible  time-varying  linear  system.  Hamon  and  Hannan  have  previously 
described  an  approximate  ML  solution  to  the  time-delay-estimation  problem 
for  dispersive  systems. 
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Although  the  details  of  our  derivation  differ  substantially  from 

those  In  reference  1,  the  estimation  criterion  Is  Identical.  We  seek  a 

function,  the  log-likelihood  function  1  rvA< 0( t ) , A*) ,  whose  value  Is  an 

Indicator  of  the  likelihood  that  a  hypothetical  delay  function,  D(t),  and 

attenuation  constant,  X,  caused  a  particular  received  vector  waveform,  R(t), 

& 

<  t  <  Tf.  The  function  0(t)]ML  and  the  constant  A]ML  jointly 
maximizing  ln/UD(t),X)  are  the  ML  estimates  of  d(t)  and  X,  respectively, 
when  R(t)  Is  the  received  vector  waveform. 

This  report  is  organized  as  follows:  The  log-likelihood  function 
lnA(D(t),X),  derived  in  section  2,  Is  shown  to  depend  upon  the  minimum  mean 
square  error  (MMSE)  estimators  of  s(t)  and  as(t  -  d ( t ) )  from  r(t) 
conditioned  on  given  attenuation  and  delay.  We  show  how  to  implement  these 
estimators  in  section  3.  In  section  4  we  obtain  the  four  alternative 
systems  for  computing  lnA<D(t)  ,X) .  In  section  5  we  show  that  the  general 
solution  to  the  problem  of  ML  estimation  of  d(t)  reduces  to  the  generalized 
crosscorrelator  receiver  In  reference  1  for  the  special  case  that  d(t)  is  a 
constant,  the  signal  process  is  stationary,  and  the  observation  interval  is 
long. 
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2.  DERIVATION  OF  THE  LOG-LIKELIHOOD  FUNCTION 


The  derivation  of  the  log-likelihood  function,  lnA(D(t),A),  Is  somewhat 
lengthy  and,  therefore,  has  been  separated  Into  three  parts.  The  first  part 
of  the  derivation,  described  in  subsection  2.1,  obtains  a  series  form  for 
1  oAC D( t )  ,£*)  using  the  generalized  Karhunen-Loeve  expansion.  The  result  is 
given  in  equations  (2-16),  (2-17),  and  (2-18).  The  second  and  third  parts, 
in  subsections  2.2  and  2.3,  show  that  the  series  (2-17)  and  (2-18)  can  be 
put  into  the  closed  forms  (2-19)  and  (2-36),  respectively,  that,  in  turn, 
depend  upon  the  noncausal  and  the  causal  MMSE  estimators  of  s(t)  and 
3*s(t  -  d ( t ) )  from  r(t)  conditioned  on  given  attenuation  and  delay. 


The  developments  in  subsections  2.1,  2.2,  and  2.3  are  basically 
extensions  of  the  material  in  references  17  (pp.  203-205  and  221-223)  and  19 
(pp.  22,  pp.  170-173)  from  scalar  to  vector  random  processes. 


2.1  SERIES  FORM 


The  problem  of  estimating  attenuation  and  time-varying  delay  can  be 
reframed  as  a  parameter-estimation  problem.  One  way  to  do  this  is  to 
represent  the  time-varying  delay,  d(t),  by  a  generalized  Fourier  series, 


d(t)  =  ^  ;  Ti  <  1  1  Tf  . 


(2-1) 


where 


B 
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d(t)^(t)dt  (2-2) 

i 

and  where  is  any  convenient  set  of  basis  functions  that  is 

complete  and  orthonormal  (CON)  over  the  Interval  [T^,Tf].  Because  d(t) 
is  nonrandom  but  unknown,  the  coefficients  d1 ,  d2,  d^,  ...  are 
nonrandom  but  unknown.  The  substitution  of  (2-1)  into  s(t  -  d(t))  yields  a 
function  that  depends  upon  the  basis  set  |4^(t)|  ,  the  vector  of  unknown 
coefficients  d  =  (d-j.dj,  ...),  and  t.  To  show  the  dependency  on  d 
explicitly,  we  will  denote  this  function  by  s(t;d);  that  is. 


s(t;d)  £  s ( t  di+i(t)).  (2-3) 

i*l 

It  follows  from  notation  (2-3)  that 


s(t;0)  «  s(t) 


(2-4) 


We  now  write  r(t)  of  (1-2)  as 

r(t)  =  s(t;d,a)  +  w(t)  ,  (2-5) 

where 

r(t)  =  ( r1 ( t)  r2( t) )T  ,  (2-6a) 

s(t;d,a)  =  ( s < t ; 0 )  as(t;d))T  , 


(2-6b) 
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w(t)  =  (w^t)  w2(t))T  .  (2-6c) 

The  problem  of  estimating  the  unknown  delay  d(t)  and  relative 

attenuation  constant  3  from  r(t)  In  (1-2)  Is  equivalent  to  that  of 

estimating  the  unknown  vector  d  and  scalar  ‘a  from  r(t)  In  (2-6).  We  have 

considered  the  generalized  Fourier  series  representation  of  d(t)  in  (2-1) 

because  It  is  both  well-known  and  general.  Other  techniques  for  representing 

d(t)  as  a  vector  may  be  preferred  in  a  particular  application.  For  example, 

2 

if  it  is  known  a  priori  that  d(t)  =  ag  +•  a^  +  a^t  for  <  t  < 

(where  a^,  a^ ,  and  a2  are  unknown),  one  can  set  d  =  (a^,  a^ 

a^)T  to  obtain  a  four-parameter-estimation  problem  involving  the  unknown 
attenuation  3  and  the  physically  meaningful  constants  aQ,  a1 ,  and  a2. 

Notice  that  if  d(t)  Is  a  known  function,  D(t),  and  if  3  is  a  known 
constant  /C,  then  as(t  -  d( t ) )  is  related  to  s(t)  simply  by  a  known  linear 
transformation,  L  ~{s(t)}  =  3s(t  -  0(t)).  Therefore,  for  given  3  and 

U  V  L /  ,A 

d(t),  the  signal  s(t)  and  its  delayed,  attenuated  version  3s(t  -  d(t))  are 
jointly  Gaussian  random  processes.  Let  D  be  the  vector  of  coefficients 
corresponding  D(t).  Then,  for  d  =  D  and  3  =  A",  r(t)  in  (2-6)  is  a  Gaussian 
random  vector  process  having  mean  zero  and  2X2  matrix  covariance  function, 


Kr;d,i(t.u.;0,A)  =  E{r(t)rT(u)|d  =  D ,a  =  a[ 

=  E|s(t;D,A)sT(u;D,A)|  +  Ejw(t)wT(u)} 

N 

=  Ks.d~(t,u;D,A)  +  I  6(t  -  u)  ,  (2-7) 


7 
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where  I  Is  the  2X2  Identity  matrix. 


We  proceed  by  representing  vector  process  r(t)  as  an  Infinite 


dimensional  vector  r  using  the  generalized  Karhunen-Loeve  expansion 


(references  17  (pp.  221-223)  and  20).  We  define 


ii 

rN(t)  ^  ^  r.0.(t;D,A)  , 


(2-8) 


where 


fTf  T  - 

r i  =J  0.  (t;0,A)r(t)dt  ;  1  -  1.  2 . 


(2-9) 


and  where  the  0. (t;D,A)  are  the  normalized  vector  eigenfunctions  of  the 
matrix  covariance  function  X  .  ~(t,u;D,X) .  We  assume  that  (0.-(t;D,A)} 

S ,Q ta  1 


is  complete.  The  normalized  vector  eigenfunctions  are  two-element-col umn 


vectors  that  satisfy  the  equations 


.  rTf  . 

S(9-A)5i(t;O.A)  =J  Ks.d~(t,u;D,A)0i(u;D,A 


A)  du  ;  T^<  t  < 


(2-10) 


rTf  t  - 

J  0i  (t;O,A)0j(t;g,A)  dt  =  6.^  . 


(2-11) 


where  x.(0,X)  is  the  (scalar)  eigenvalue  associated  with  0.(t;D,A). 


With  the  0.(t;D,X)  so  specified,  it  follows  that 


8 
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r(t)  «  1.1. m.  r^t)  .  (2-12) 

N  -»  oo 

This  Is  the  generalized  Karhunen-loeve  expansion  of  r(t). 

Because  {^(tiO.A)}  Is  complete,  one  can  represent  r(t)  by  (2-12) 
(using  (2-8)  and  (2-9))  for  hypothetical  or  assumed  values  of  D  and  A.  It 
Is  easy  to  show  that  If  the  assumed  values  of  D  and  A  are  the  true  values  of 
the  unknown  quantities  d  and  S',  respectively,  then  the  n’s  in  (2-9)  are 
statistically  Independent  Gaussian  random  variables  having  zero-  means  and 
variances. 


E<{ri 2|d  =  0,  a  *  A|  =  \.(0.A)  ♦  ^  ;  1  -  1,  2 . N  .  (2-13) 

The  joint  probability  density  function  of  the  r^  conditioned  on  d  -  D  and 
a  *  A  is,  therefore. 


prN;d,a(-N:-,A) 


N 

7T 


i*i  57 

^2*[\.(0,A)  +  f] 


exp 


2[Ki(0,A)  +  f] 


(2-14) 


where  =  (r-|  r2  •••  rfj)T  and  =  (R1  To 

obtain  the  likelihood  function  associated  with  r(t),  we  take  the  logarithm 
and  the  limit  N  -» «  .  This  leads  to  a  convergence  difficulty  that  can  be 
bypassed  in  the  usual  way  (see  reference  17,  p.  274)  of  dividing  (2-14)  by 
the  function 
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N  1 

7T  —  exp 
i«l 


f(RN)  =  7 T  —  exp  -  jr— 


(2-15) 


The  result  of  this  division  is  still  a  likelihood  function  because  f(R.N) 
does  not  depend  upon  0  or  /C.  After  dividing  (2-14)  by  (2-15),  we  take  the 
logarithm  and  the  limit  N  The  result  is  the  log-likelihood  function 


lnA(D,A)  =  iR( D.£)  +  iB(0.X)  . 


(2-16) 


where 


•  / n  7 \  a  J.  V  ' 

(  Q  *  A )  -  u  /  _ 

0  W]  (O.A)  +  N 


Xi  (^'A)  R2 


(2-17) 


V0-*5  “  ~  \  1n[1  +  N^  y0’A>) 

i=l  0 


(2-18) 


2.2  CLOSED  FORM  FOR  DATA  DEPENDENT  TERM  lR(D,A) 


The  first  term  in  (2-16),  iR(0,X),  can  be  written  as 


yo.A)  = 


I  r  I, 

iff 


R  (t)H  (t.v;0,A)R(v)dtdv  , 


(2-19) 


V 


TR  7665 


Hn(t,v;O.A)  £ 


y;  X^D.A)  g^tiD.Ax/fviP.A)  ; 
1-1  X^O.A)  +  Nq/2 


T.  <  t.v  <  T.  , 


(2-20) 


and  R(t)  is  the  sample  vector  function  of  r(t)  corresponding  to  the  vector  R: 


oo 

R(t)  =  y  Ri«i(t;D,A)  , 


(2-21) 


rTf 

>1  =J  li^tjO.A 


A)R(t)dt  . 


(2-22) 


Equation  (2-19)  can  be  verified  by  substituting  H  (t,v;D,A)  in  (2.20) 


into  (2-19)  and  using  (2-22).  An  interpretation  of  (2-19)  is  obtained  by 


considering  the  following  noncausal  linear  estimate  of  s(t;D,A): 


L(t;0,A)  = 


fTf 

J  Hn(t,v;0,A)r(v)dv,  T.  <  t  <  Tf  , 


(2-23) 


where  the  subscript  "n"  is  used  to  denote  a  noncausal  estimate  or  system. 


It  follows  from  (2-23)  and  (2-7)  that 


*  •  •  vv.v . v.v.v.v v.-cv^wv  •'*•**•••• 
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E{[s(t;d,a)  -  in(t;0,A]rT(u)|d  =  D.a  =  A^=  Ks;d~(t,u;0,A) 


N 

-  y  «n(t.u;D.A) 


(2-24) 

20 

According  to  the  matrix  version  of  Mercer's  Theorem: 


T 

J 


Hn(t,v;0,A)Ks.d~(v,u;0,A)dv 


Ti  <- 


t,u  <  T 


f  * 


Yj  X.(D,A)0.(t;D,A)e.T(u;O,A) 
i=l 


Ti  <  *• 


<  Tf 


(2-25) 


By  substituting  (2-25)  and  (2-20)  into  the  right-hand  side  of  (2-24)  and 
using  (2-11)  one  obtains 


KS;dfj(t.u;0.A)  -  f  H 


n(t.u;0,A)  -I  Hn(t,v;D,A)K  d>j(vtu;D,A)dv 


T. 
' l 


0  ;  T.  <  t.,  u  <  Tf  . 


( 2 -2  6a ) 


which  yields 


£{[s(t;d,a)  -  fn(t;0,A) ]rT(u)  J  d  =  0,a  =  a|=  0  ;  T.  <  t.u  <  Tf.  (2-266) 


Therefore,  if  D  and  A  are  the  true  values  of  d  and  a-,  then  the  estimation 
error 


en(t;d,a.0,A)  =  s(t;d,a)  -  fn(t;D,A) 


(2-27) 
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Is  orthogonal  to  r(u),  <  u  <  Tf.  Consequently  (see  reference  21, 

p.  390)  i^tjO.A)  In  (2-23)  Is  the  noncausal  point  LMMSE  estimate  of 
s.(t;d,T)  from  r(u),  <  u  <  Tf,  given  that  d  =  D  and  a"  =  A  and 

H^t.vjO.A)  is  the  Impulse  response  of  the  noncausal  point  LMMSE 
estimator.  Note  that  the  2-by-2  matrix  function  H^t.vjD.A)  is  the 

solution  to  equation  (2-26a).  As  will  be  described  in  section  4,  the 
substitution  of  (2-23)  into  (2-19)  results  in  a  vector  estimator-correlator 
realization  for  lR(0,/f). 

The  error  convariance  matrix  of  the  noncausal  point  LMMSE  estimate  of 
s(t;d,a)  conditioned  on  d  =  0  and  H  =  A  is 

E„(t;0.»)  S  E{e„(t;d.S.D.A)enT(t:d.5.0.A)|d  *  0,  l  =  a}  .  (3-28) 

By  substituting  (2-27)  into  (2-28)  and  using  (2-23)  and  (2-26),  one  obtains 

N 

En(t;D.A)  =  ^  Hn(t,t;D,A)  .  (2-29) 

2.3  CLOSED  FORM  FOR  BIAS  TERM  i0(D,A) 

D 

The  term  IgtD.A-)  in  (2-18)  can  be  written  in  closed  form  by  noting 
that,  according  to  (2-10)  and  (2-11),  the  eigenvalues  \.(D,A)  and  the 
vector  eigenfunctions  0^(t;0,A)  depend  upon  the  final  observation  time 
Tf.  To  indicate  this  dependency,  we  write 


M(D.A)  =  Xi(0,A,Tf) 


(2-30) 
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ii(t;D,A)  =  ti(t;D,A.Tf)  . 


It  follows  that  !o(0,A)  In  (2.18)  can  be  rewritten  as 


(2-31) 


lB(0.A)  = 


-  \ J  dt  ^  £  TnCI  +  J-  X.(D,A,i 


.  f1  f  «  [d\. (0,A,t) ]/dt 

r  dt  /  - - - :: —  • 

o-'  frf  1  +  (2/N  )Xi(D,A,t) 


(2-32) 


where  x.(0,A,T.)  =  0.  It  can  be  shown  by  a  straightf orward  extension  of 
the  derivation  in  reference  17,  pp.  204-205,  that 


d\.(0,A,t) 


=  Xi(O,A,t)Tr|0i(t;O,A,t)e.T(t;D,A,t)j>  . 


(2-33) 


When  we  use  (2-33)  and  the  fact  that  Tr  {•}  is  a  linear  operator,  (2-32) 


becomes 


yo.»)  -  -  ? 


('  J±  - 

J  i-l  v,(  O.A.t)  t  N  n 
i 


0.(t;O,A,t)0.  (t;D,A,t)|dt  . 


(2-34) 


A  closed  form  for  the  quantity  in  the  braces  in  (2-34)  is  recognized  by 
rewriting  (2-20)  with  the  notation  of  (2-30)  and  (2-31): 
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X .  ( 0 ,  A , T, )  t 

Hn(t.v;D.X.Tf)  "  2J— 4  —  li(t:M.Tf)iiT(v;D.y.Tf)  . 


(2-35) 


The  quantity  in  the  braces  in  (2-34)  is  Hn(t,t;0,A,t) .  Since 
H  (t,v;D,A,t)  is  the  matrix  impulse  response  of  the  point  LMMSE  estimator 


of  s(t;d,lT) 

from 

r(v).  T.  < 

4-* 

VI 

> 

given 

d 

=  D 

and  K 

=  a. 

then 

^(t.viO.A.t) 

is. 

by  definition,  the 

matrix 

impulse 

response  of 

the 

causal  point 

LMMSE 

estimator  of 

s(t;d  ,t) 

from  r( v) , 

given 

d  =  D 

and  a" 

=  A. 

If  we  denote 

the  causal  matrix 

impulse 

response 

by 

H  (t. 

v;D,A) , 

then 

(2-34) 

becomes 


1  f  f 

lB(D,A)  =  -  Tr(Hc(t,t;D,A)]  dt  .  (2-36) 

T. 

All  the  previous  equations  describing  noncausal  estimation  of 
s(t;d,a*)  from  r(v)  describe  causal  estimation  of  s(t;d,aO  from  r(v)  when  t 
Is  substituted  for  T^.  In  particular,  with  =  t,  equation  (2-29) 
describes  the  error  covariance  matrix  of  the  causal  LMMSE  estimate  of 
s ( t ; d , T) .  given  d  =  D  and  "a"  =  AT: 


Ec(t ;D.A)  =  f  Hc(t,t;0,A)  . 


(2-37) 


The  substitution  of  (2-37)  into  (2-36)  yields  an  alternative  expression 


for  lB(0,A): 


yq.A)  ■  -  jjl 
0 


fTf 

J  Tr(Ec(t;0,A 


A)]  dt  . 


(2-38) 


T . 
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3.  THE  MATRIX  IMPULSE  RESPONSE  H  (t,v;0,I) 

— n  — 

In  this  section  we  derive  a  simple  explicit  form  for  the  matrix  impulse 
response  H^t, v;D,X) .  It  is  relatively  difficult  to  obtain  this  form  by 
solving  equation  (2-26a).  The  constructive  approach  of  subsection  3.2  has 
the  advantage  of  being  both  mathematically  and  conceptually  simple.  Before 
proceeding  with  the  constructive  solution,  it  will  be  helpful  to  derive  the 
explicit  form  for  the  inverse  of  the  operator  (1-4). 

3.1  INVERSE  OF  OPERATOR  (1-4) 

By  definition,  the  inverse  operator  satisfies 


|s ( t ; 0 , A) }  =  LD( t) , A  ^S(t  '  0(t))f  =  S(t)  ■  (3-1) 

Let  v(t)  be  an  arbitrary  waveform  and  try  an  inverse  having  the  form 

LD(t).A  lv(t)l  =  ~  v(B(t))  .  (3-2) 

where  I3(t)  is  to  be  determined.  Define 

f(t)  =  t  -  D(t)  (3-3) 

so  that  by  the  definition  (3-1) 


(3-4a) 


L0(t)  ,A  l“<f<t>»l  • 
Replacing  t  by  f(t)  in  (3-2)  gives 


L0(t),Aiv(f(t))i=  ~  v(S(f(t)))  *  (3_4b) 

which,  when  compared  with  (3-4a),  yields 

B(f (t))  =  t  .  (3-5) 

Therefore,  the  inverse  operator  is  given  by  (3-2),  where  B(«)  is  the  inverse 


of 

the  function 

f(t). 

S1nce  L0(t).Jf-> 

is 

invertible,  the  function 

f(t) 

is 

one-to-one. 

We 

now  proceed  to 

the 

constructive  derivation 

of 

yt.vjD.ff). 

3.2  CONSTRUCTIVE  DERIVATION  OF  ^(t.vjD.A) 

The  first  step  in  the  derivation  of  Hn(t,v;D,A)  is  to  ( noncausal ly) 
transform  r(t),  T.  <  t  <  T  ,  into  the  vector  process  r'(u), 
f ( T . )  <  u  <  T^ ,  where 

r~-i 

A  r  (B(u) ) 

r'(u)  = 

0 

r^u)  + 

r  -j  ( ci )  - 


;  f(T.)  <  u  <  T1  ,  ( 3-6a) 

A-1 r2(0( u) ) 

;  T.  <  u  <  f(Tf)  ,  (3-6b) 

A_1  r2((3(u) ) 


TR  7665 


;  f(Tf)  <  u  <  Tf  ,  (3-6C) 

and  f(u)  is  defined  in  (3-3).  In  (3-6),  X  can  be  regarded  as  an  assumed 

value  for  the  unknown  relative  attenuation  constant  af,  and  D(t)  as  an 

assumed  function  for  the  unknown  delay  function  d(t).  (We  naturally 

require  D(t)  >  0,  which  implies  that  f(t)  <  t.)  Notice  that  (3-6b)  assumes 

that  f(Tf)  is  not  less  than  T..  This  is  equivalent  to  the  assumption 

that  the  signal  delay  does  not  exceed  the  observation  interval.  Since  this 

assumption  is  likely  to  be  met  in  most  applications,  we  will  retain  it  in 

the  following.  It  is  not  hard  to  generalize  our  results  to  include  the  case 

of  f (T  )  <  T.. 
f  i 

The  transformation  r(t)  -»  r'(u)  is  illustrated  in  figure  3-1,  where, 
for  simplicity  in  interpretation,  the  noise  processes  w^t)  and  w2(t) 

have  been  drawn  as  small  ripples.  A  system  block  diagram  for  the 
transformation  is  shown  in  figure  3-2.  An  examination  of  equation  (3-6), 
figure  3-1,  or  figure  3-2  will  reveal  that  the  transformation  from  r(t)  to 
r'(u)  is  linear  and  invertible.  Thus,  r(t),  T..  <  t  <  Tf,  can  be 

recovered  from  r'(u),  f(T^)  <  u  <  T^,  using  a  linear  transformation.  It 
follows  from  the  reversibility  theorem  (reference  17,  p.  289)  that  the 
noncausal  IMMSE  estimate  ^(t^.A)  in  equation  (2-23),  given  d  =  D  and 

a"  =  A,  can  be  obtained  from  r'(u).  Before  describing  the  structure  of  the 
IMMSE  estimator,  it  will  be  helpful  to  observe  that  if  d  =  D  and  a '  -  1\, 

then,  from  equations  (2-5),  (2-6),  and  (3-6), 


r'(u) 


^(u) 
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r'(u)  =  0  ;  u  <  f(T-j), 


S(U) 

r'(u)  *  + 

0  . 


["’'“’I  ; 

ln2(u)J 


f(Ti)  <  u  <  Tf  , 


(3-7a) 


( 3— 7b) 


r'(u)  =  0  ;  Tf  <  u  , 


( 3— 7c ) 


#here 


n-j  ( u ) 


.n2(U). 


A  w2(0(u)) 


;  f^)  <  u  <  Ti, 


( 3-8a) 


Vu) 

_n2(u)_ 


W-|  ( u )  +  A  w2(0(u)) 

\  -1  :  Ti  <  u  i  f(Tf> 

2  w^u)  -  A  'w2(S(u)) 


( 3  —8b ) 


"nl  (u) 


n  2  ( u ) 


W^u) 


;  f(Tf)  <  U  <  Tf  . 


(3-8c) 


We  present  the  form  of  the  LMMSE  estimator  of  s(t),  f(T.)  <  t  <  Tf. 


in  the  following  theorem: 


Theorem.  The  noncausal  point  LMMSE  estimator  of  s(t)  from  r'(u), 
f(T..)  i  t,  u  <  Tf,  conditioned  on  d  =  D  and  a"  =  X,  is  given  by  the 
system  in  figure  3-3,  where  f(t,u;0,X)  is  the  impulse  response  of  the 
noncausal  point  LMMSE  estimator  n,  (t)  of  n,(t)  from  n0(u). 
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±'(U) 


r-,'  (u) 


r4(u) 


gn(t,u;D,A) 


f(t,u;D,A) 


Figure  3-3.  Structure  of  the  Noncausal  Conditional  IMMSE  Estimator  of 
s(t)  From  r 1  (u) ,  f (T^ )  <  u  <  Tf  (When  d  =  D  and  a  =  S',  then  x(t) 
equals  the  noncausaj,  LMMSE  estimate^of  s(t).  The  Impulse 
responses  f(t,u;0,A)  and  gn(t,u;D,A)  are  defined  by  the 
theorem  in  section  3.  f(t,u;0,7C)  is  specified  by 

equations  (3-30)  and  (3-31),  and  gn(t,u;D,A)  is 
specified  by  equations  (3-44)  and  (3-41). 
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^(t)  *  J  f(t,u;0,A)n2(u)du  ;  f (Ti )  <  t  <  Tf  .  (3-9) 

f(f,) 

and  gn(t,u;0,S')  Is  the  Impulse  response  of  the  noncausal  point  LMMSE 
estimator  ^n(t)  of  s(t)  from  s(u)  +  n^(u)  -  ^(u), 


^n(t)  *  J^g^t.ujD.A)  [s(u)  +  n^u)  -  ^(u)]du  ;  f(T.)  <  t  <  Tf 

fdi) 


(3-10) 


Proof.  According  to  the  orthogonality  principle  (reference  21,  p.  390) 
a  linear  functional  p  =  L[g]  is  the  LMMSE  estimate  of  a  random  variable  p 
from  data  vector  g(§)  §cD  (where  0  is  the  domain  of  the  data)  if  and  only  if 
the  estimation  error  p  -  p  is  orthogonal  to  g  for  all  $cD, 


E  j  (P  -  P)g($)}  =  0  ;  |C0  .  (3-11) 

Therefore,  a  necessary  and  sufficient  condition  that  ^(t)  be  the  LMMSE 
estimate  of  s(t)  from  r'(u),  given  that  d  =  D  and  If  =  X,  is  that  the  vector 


v(t.u)  =  E  | [s(t)  -  £n(t)]  r’(u)  |  d  =  0,  I  =  A  |  (3-12) 

be  identically  zero  for  f(T.)  <  t,u  <  T^.  Using  (3-7),  we  note  that  the 

components  of  y  (t,u)  are 

v^t.u)  =  Ej[s(t)  -  ^n(t)]  [s(u)  +  n^u)]}  (3-13) 

and 
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v2(t,u)  =  E  \  [S(t)  -  $n(t)]  n2(u)| 


for  f (T1 )  <  t.u  <  Tf 


By  the  definitions  of  f(t,u;D,A)  and  gn(t,u;D,A), 


EjCn^t)  -  fr1(t)]n2(u)}  =  0  ;  f(T.)  <  t.u  <  Tf 


(3-14) 


(3-15) 


Ej[s(t)  -  ^n(t)][s(u)  +  n] (u)  -  ^(u)]|  =  0  ;  f(T.)  <  t.u  <  Tf. 


(3.16) 


Recall  that  the  signal  process  s(t)  is  orthogonal  to  the  white  noise 
processes  w^t)  and  w2(t).  Therefore,  s(t)  is  also  orthogonal  to  the 
noise  processes  n^t)  and  n?(t)  defined  in  (3-8). 

It  follows  that  (3-14)  simplies  to 


v?(t.u)  =  -E{$n(t)n2(u)} 


which,  with  the  aid  of  (3-10),  becomes 


( 3-1 7a) 


v2(t,u)  =  -Ej  J  gn(t,o;0,A)  [s(o)  *  n }{a)  -  ^ (a)  ]do]n2(u) 

hTi) 

rTf 

=  J  gn(t,o;0,A)E{[n1(o)  -  ft1(o)]n2(u)}  do 
fdi) 


=  0  ;  f ( T . )  <  t,  u  <  Tt  . 


(3-1 7b) 
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where  the  last  step  follows  from  (3-15).  The  result  that  v^t.u)  in 
(3-14)  equals  zero  and  the  fact  that  ^(t)  depends  linearly  upon  n^(t) 
together  imply 

E{[s(t)  -  ^n(t)]^(u)}  =  0  ;  fd^  <  t.  u  <  Tf  .  (3-18) 

Substracting  (3-18)  from  (3-13)  leads  to 

Vl(t.u)  =  E{[s(t)  -  ^n(t)]  [s(u)  +  n^u)  -  ^(u)]}.  (3-19) 

Comparing  (3-19)  with  (3-16),  we  see  that 
v^t.u)  =  0  ;  f(T.)  <  t.  u  <  Tf  .  (3-20) 

This  completes  the  proof. 

The  LMMSE  estimator  of  a"s(t  -  d(t)),  T^  <  t  <  Tf,  from  r(u), 

T^  <  u  <  Tp  conditioned  on  d  =  D  and  =  ft,  follows  easily  from  the 
fact  that  fts(t  -  d(t))  is  a  linear  transf ormation  of  s(t).  Because  all 
available  data  have  been  used  to  obtain  £  (t),  f(T.)  <  t  <  T.,  the 
noncausal  LMMSE  estimate  of  fts(t  -  d(t)),  given  d  =  D  and  1  -  ft,  is  simply 
the  scaled  and  delayed  version  of  ^n(t)  in  (3-10),  namely,  ft£  (t  -  D(t)). 

The  explicit  form  for  f(t,u;D,ft)  follows  easily  from  (3-9)  and  (3-15), 


which  together  imply 
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fJf 

E{n1(t)n2(u)J  =  J  f(t,0;O,'X)E{n2(o)n2(u)  J  do  ,  f (T^ )  <  t,  u  <  Tf  . 


f(T.) 


This  can  be  simplified  by  using  (3-8)  and  (1-5),  which  imply 


(3-21) 


E{n2(t)n2(u) 


N,  r  ?  1 

j  |i(t  -  u)  +  <S(B(t)  -  B(u))J ;  T.  <  u  <  f(Tf) 


0  ;  otherwise 


(3-22) 


E{n1(t)n2(u) J 


N  . 

[i(t  -  u)  -  A  6(B(t)  -  B(u))J;  T.  <  u  <  f(Tf)  . 


10;  otherwise 


(3-23) 


Since  B(t)  is  a  one-to-one  function,  then 


4(6(t)  -  fl(u) )  =  — - -  4(t  -  u)  , 

|fl(u)| 


(3-24) 


where  the  dot  denotes  the  derivative  of  a  function.  It  follows  from  (3-3) 
and  (3-5)  that 


0(u)  = 


1  -  D(B(u)) 


(3-25) 


Note  that  s(t  -  0( t ) )  is  locally  reversed  in  time  where  D(t)  >  1  and 
frozen  in  time  where  D(t)  =  1.  Therefore,  it  is  reasonable  to  define  D(t) 
as  a  valid  delay  function  if  and  only  if 


D( t)  <  1  , 


(3-26) 


s'  \VV. •'.V.V.V.  V  \  .-.  '  ••  . 
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where  the  equality  holds  only  at  Isolated  values  of  t.  It  is  easy  to  see 
that  the  above  condition  guarantees  that  f(t)  is  invertible.  By  combining 
(3-24),  (3-25),  and  (3-26),  we  obtain 


«(0(t)  -  B(u) )  *  [1  -  D(0(u) ) ]  4(t  -  u)  . 


Therefore,  (3-22)  and  (3-23)  become,  respectively. 


(3-27) 


Ejn2(t)n2(u)|  = 


_o  -|  +  l  -  D(B(uli  4(t  -  u)  ;  T.  <  t  <  f(Tf) 
A  c 


0  ;  otherwise 


(3-28) 


s*  1  - 


E|n1(t)n2(u)}  = 


i(t-u)  ;  T.  <  u  £  f ( Tf ) 
A 


0  ;  otherwise 


(3-29) 


The  substitution  of  (3-28)  and  (3-29)  into  (3-21)  yields 


f(t,u;0,£)  =  k(u)  5(t  -  u)  , 


(3-30) 


where 


^ - U - ^LLLL  ;  T  <  u  <  f(T  ) 

A^  ♦  [1  -  D(B(u) ) ]  1 


0  ;  otherwise. 


(3-31 


Therefore , 
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^(t)  =  k(t)n2(t)  . 


(3-32) 


Note  that  the  MMSE  estimate  of  n^t)  from  n^ft)  requires  only 
multiplication  of  n^(t)  by  a  time-varying  gain. 


It  is  interesting  to  observe  from  (3-29)  that  if 


0(t)  =  Dq  +  (1  -  TT)t 


then  n^(t)  and  n ^ ( t )  become  statistically  independent  and 


(3-33) 


^(t)  =  0  . 


(3-34) 


Equation  (3-33)  is  a  necessary  and  sufficient  condition  for  (3-34).  A 
sufficient  condition  arises  when  the  delay  is  constant, 


0(t)  =  Dq 


and  the  magnitude  of  the  attenuation  constant  is  unity. 


A  =  +1  . 


( 3-35a ) 


( 3 —35b ) 


These  results  are  a  consequence  of  the  fact  that  the  statistics  of  w^(t) 
are  unchanged  by  the  inverse  operator  (3-1)  when  D ( t )  and  A-  satisfy  (3-34). 
The  point  is  that  if  d(t)  and  7  are  known  a  priori  to  satisfy  (3-33)  then 
the  hypothetical  quantities  0 ( t )  and  "a  can  also  be  assumed  to  satisfy 
(3-3).  This  results  in  a  simplified  receiver  because  under  these  conditions 
k(u)  is  identically  zero. 
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An  equation  specifying  gn(t,u;D,A)  can  be  obtained  by  using  the  fact 
that  gn(t,u;0,/C)  is  the  LMMSE  estimator  of  s(t),  f(T^)  <  t  <  Tp  from 

z(u)  =  s(u)  +  n(u)  ;  f (Ti )  <  u  <  Tf  ,  (3-36) 


where 


n(u)  =  n-|(u)  -  fy(u)  ;  f  ( Ti )  <  u  <  Tf  .  (3-37) 

The  noise  process  n(u)  Is  zero  mean  and  uncorrelated  with  the  signal  process 
s(u).  Its  covariance  function  is 

E|n(t)n(u)f  =  E {[n] (t)  -  ^  (t)  (u)  -  fyu)]} 

-  Ejtn^t)  -  n^  (t)  ]n]  (u)j 

*  Ejn^tjn^u)}-  k(t)E{n2(t)n1(u)|.  (3-38) 

By  a  derivation  similar  to  that  leading  to  (3-28)  one  finds 


E|n] (t)n] (u)|  = 


N 

[1  -  6(fl(u))]  4(t  -  u)  ;  f  ( T  - )  <  t  <  T 
2A^  i-i 


1  + 


1  -  D(B( u) ) 


A2 


4(t  -  u)  ;  T.  <  t  <  f(Tf; 


N 

-f  4(t  -  u)  ;  f(Tf)  <  t  <  Tf  . 


Combining  (3-38),  (3-39),  and  (3-31)  gives 


(3-39) 


E  }n( t) n( u) }  =  0(u)4(t  -  u)  ;  f(T.)  <  u  <  T 


(3-40) 
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where 


Q(u) 


N 

~Y  n  -  6(B(u))]  ;  fd^  <  u  <  T1 


No  [1  -  D(P(u),J -  ;  T.  <  u  <  f (Tf ) 

2[r  +  [1  -  D(B(u))]j 


[-|  ;  f(Tf)  <  U  <  Tf  , 


(3-41) 


and  it  follows  that  n(t)  is  nonstationary  white  noise.  The  equation 
specifying  gn(t,u;0,/f)  is  now  obtained  by  substituting 


£„<t) 


/' 


9n( t,o;D,A) z(o)do 


fd^  <  t  <  Tf 


into  the  orthogonal ity  condition  , 


(3-42) 


E|[s(t)  -  's  n  ( t )  ]  z(u){  =  0  ;  f(T.)  <  t.u  <  Tf.  (3-43) 

This  leads  directly  to 


Rs(t.u)  =  J  9n(t,<j;D,A)  Rs(o,u)  do 
f<Ti) 

*  0(u)gn(t,u;0.A)  ;  f(T.)  <  t.u  <  Tf  .  (3-44) 

Note  that  the  problem  of  finding  gn(t,u;D,X)  is  equivalent  to  the 
problem  of  deriving  the  noncausal  LMMSE  estimator  of  s(t)  from  s  ( t )  «-  n(t), 
where  the  process  n ( t )  is  nonstationary  white  noise  uncorrelated  with  s(t). 
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If  s(t)  is  a  state  representable  process,  then  ?n(t)  can  be  obtained  using 
optimal  linear  smoothing  (reference  22,  chapter  5). 


3.3  EXPLICIT  FORM  FOR  ENTRIES  IN  H  (t,u;D,A) 


We  have  now  specified  the  structure  of  H^t.ujD.A).  This  structure 

Is  shown  in  figure  3-4,  where  the  filter  gn(t,u;D,/T)  is  given  by  the 

solution  to  (3-44)  and  where  k(t)  is  given  by  (3-31).  Using  this  structure, 

we  now  derive  the  explicit  form  for  the  individual  entries  h. j(t,u;D,/T)  in 

H  (t,u;0,A).  This  can  be  done  be  noting  that,  by  definition,  the  output 
— n  — 

of  H?)(t,u;D,/T)  is  (see  figure  3-4) 


x(t)  =  h11(t,u;0,S)r1(u)du  +  h12(t,u;0,S)r2(u)du  (3-45) 


y(t)  =  I  h21(t 


rTf 

,u;0,A)r^(u)du  +  I 


h22(t,u;D,A)r2(u)  du  .  (3-46) 


On  the  other  hand,  by  tracing  the  signals  through  the  system  in  figure  3-4, 

we  can  express  x(t)  and  y(t)  in  terms  of  g  ( t ,  u ;  D .  A*)  and  k(u).  To  keep 

n  - 

the  notation  simple,  we  will  subsequently  write  g  (t,u;D,A)  as  gn(t,u). 

An  examination  of  figure  3-4  with  the  aid  of  figure  3-2  and  equation  (3-6) 


yields,  after  a  little  labor. 
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/  *  *n(t-u>  j  r; 


(S(u))du 


J  9n(t.u)  j  [1  -  k(u)]  r1(u)du 


rf(Tf) 

gn<t,u)  h  l>  *  k(u)]  r  (0(u))du 
J-  "  2A  4 


gn(t,u)r(u)du1  . 


f  (Tf ) 


(3-47) 


By  changing  variables  in  the  first  and  third  integrals  (set  a  =  B(u)) 
(3-47)  becomes,  after  a  little  more  labor. 


(*0(Ti} 

=  J  V t.o  - 


x(t)  = 


0(o))  -  r_(o)  [1  -  D(«»)]do 
A 


rf(Tf} 

+  J  2 


k(u)]r1(u)du 


*  ('  v*. 


B(T.) 


o  -  0(o))  —  [1  +  k(o  -  D(o))]r  (o)  [1  -  0(o) ]do 
2A  £ 


♦  (' 


g„(t,u)r.(u)du  . 


(3-48) 


..jr’J'Jy 
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By  comparing  (3-48)  with  (3-45)  we  obtain  for  <  t,  u  <  Tf 


h-|  i  ( t ,  u  ;0,  A )  = 


gn(t-u)2C1  "  k(u)]  ;  Ti  -  u  -  f(Tf) 


9n(t.u)  ;  f(Tf)  <  u  <  Tf 


(3-49) 


h12(t,u;D7A)  = 


gn(t,u  -  0(u) )~[1  -  D(u) ]  ;  T.  <  u  <  B(T.) 
A 


g_(t,u  -  D( u) )  — [1  +  k(u  -  D(u))][l  -  0(a)]  ;  B(T.)  <  u  <  Tf 
l  n  2A  1 


(3-50) 


We  note  that  since  a  -  D(o)  =  f(o)  then  from  (3-5)  and  (3-31) 


k(«  -  0(0))  =  %  - 

«  +  [1  -  0(a) ] 


(3-51) 


The  formulas  for  h21 (t,u;0,A)  and  h22(t,u;D,A)  in  (3-46)  can  be  obtained 
easily  by  noting  from  figure  3-4  that 


y(t)  -  Ax(t  -  D(t) )  , 


(3-52) 


which  leads  to 


h21(t.u;D,A)  =  Ahn(t  -  D( t) , u ; 0 , A) 


h22(t,u;O.X)  =  Xh12(t  -  D(t),u;0,A)  , 


where  T  <  t.u  <  T,. 


(3-53) 


(3-54) 
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3.4  EXPLICIT  FORM  FOR  BIAS  !_  (D,A) 

D 


We  can  obtain  the  explicit  form  for  the  entries  of  the  causal  matrix 

Impulse  response  Hc(t,u;D,A)  by  noting  that  if  Tf  =  t  then  none  of  the 

data  r(v),  <  v  <  Jft  are  future  data.  Thus,  for  Tf  =  t  and  for  d  = 

0  and  aT  =  X,  H^t.vjO./C)  becomes  the  impulse  response  Hc(t,v;D,X)  of  the 

causal  LMMSE  estimate  ic(t;0,/T)  of  s(t;d,a)  from  r(v),  T.  <  v  <  t,  and 

g  (t.o)  In  (3-42),  (3-43),  and  (3-44)  becomes  the  causal  LMMSE  estimator 
n 

g  (t,«j)  of  s(t)  from  z(o),  f(T^)  <  a  <  t.  Therefore,  we  can  obtain  the 
components  of  H  (t, v;0,7f)  by  replacing  gn(t,u)  in  equations  (3-49), 
(3-50),  (3-53),  and  (3-54)  with  g  (t,u),  where  gc(t,u)  is  the  solution 
to  (3-44)  for  Tf  =  t,  with  gc(t,u)  =  0  for  t  <  u.  With  H.c(t,v;D,A)  so 
determined,  lod),^)  and  fL(t;D,/T)  can  be  obtained  directly  from  (2-36) 
and  (2-37),  respectively.  For  example,  by  straightforward  substitution, 
(2-36)  becomes 


lB(D.A)  =  - 


9c^a,<^2^  ~  k(°)^da 


9c(o.o)do 


f(Tf) 

f8<V 

J  X«c 


(a  -  0(a)  ,  a  -  D(  o) )  1  -  D(a)]da 

A 


■if' 

B(T . ) 


(a  -  0(o),  a  -  D(o))— [1  +  k(o  -  D(o))](l  -  D(o))]do 

2A 


(3-55) 


,t  ut  Lt’kt  |>|A 
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This  result  can  be  put  Into  a  simpler  form  by  changing  the  variable  of 
integration  in  the  last  two  integrals.  With  a  -  0(a)  ■*  a,  (3-55)  becomes 


!b(D,A) 


-2  f'  V’- 


o)do  . 


(3-56) 


The  minimum  mean  square  error  associated  with  g  (t,u)  is 


Soc(t)  =  E{(s(t)  -  £c(t))2| 


=  Rs(t,t)  - 


gc(t,u)Rs(t,u)du 


=  Q(t)gc(t.t)  . 


(3-57) 


where  the  last  step  follows  from  (3-44),  with  Tf  =  t.  An  alternative 
expression  for  I  (0,A)  can  be  obtained  by  substituting  (3-57)  into 

D  “ 

(3-56).  This  observation  is  important  because  if  s(t)  is  a  state 
representable  process,  then  |0C(T)  can  be  obtained  from  the  matrix 


Riccatti  equation  (reference  22,  chapter  4.3). 
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4.  CANONICAL  REALIZATIONS 


Our  formulation  of  the  delay  estimation  problems  leads  naturally  to 
four  canonical  realizations,  which  are  based  upon  well-known  receiver 

structures  for  the  detection  of  Gaussian  signals  in  white  Gaussian  noise 

(reference  19,  section  2.1).  Here  we  simply  point  out  the  potential 

application  of  these  structures  In  delay  estimation.  A  more  detailed 

development  and  comparison  of  these  structures,  with  the  view  to  obtaining 
practical  estimation  algorithms,  appears  to  be  a  fertile  area  for  future 
research. 


The  substitution  of  (2-23)  into  (2-19)  (with  r( . )  replaced  by  R( . ) ) 


yields 


i^D.A) 


n  J 


RT(t)fn(t;D,A)dt 


(4-1) 


The  resulting  ML  estimator  of  d(t)  and  3"  is  shown  in  figure  4-1,  where, 
following  the  terminology  of  Van  Trees,  it  is  referred  to  as  Canonical 
Realization  No.  1.  Observe  that  this  realization  is  a  vector  estimator- 
correlator  analogous  to  the  scalar  estimator-correlator  in  figure  2-2  of 
reference  19.  Here's  how  it  works:  The  system  tentatively  hypothesizes 
that  the  unknown  delay  d(t)  is  0 ( t )  and  that  the  unknown  attenuation  'S  is  A, 
where  0(t)  is  a  possible  delay  function  and  £  is  a  possible  relative 
attenuation  constant.  The  received  vector  waveform  R(t)  is  input  to  the 
noncausal  conditional  LMMSE  estimator  of  sCtid.'J),  which  is  designed  with 
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the  assumption  that  0  and  A  represent  the  true  values  of  delay  vector  d  and 

attenuation  scalar  a\  The  output  of  the  estimator  is  ^(tjIl.A)  in 

(2-23).  A  possible  realization  of  the  estimator  was  shown  in  figure  3-4. 

The  vector  correlator  then  yields  1^(0, in  (4-1),  which,  when  added  to 

1d(D,A)  of  (2-36),  (2-38).  or  (3-55),  yields  the  value  of  the 

log-likelihood  function  lnA(0(t),A)  for  the  assumed  D(t)  and  A.  This 

process  is  repeated  for  all  choices  of  D(t)  and  A  that  are  possible  for  the 

application  in  question.  The  particular  D(t)  and  a  jointly  maximizing 

a  a 

lnA(D(t),A)  are  the  ML  estimates  0(t)]ML  and  ft]HL  of  d(t)  and  a. 

An  alternative  form  for  Canonical  Realization  No.  1  can  be  obtained  by 
noting  that  the  lower  integrator  output  in  figure  4-1  can  be  written  as 

i  rTf- 

y(t)  r2(t)  dt  =  A  x ( t  -  D(t) )  r2(t)  dt 

A  x ( o )  r2(B(o))  [1  -  D(B(o) ) ]  da  . 

(4-2) 

The  process  A-1  r? ( B( u ) )  is  available  at  the  output  of  the  inverse 
operator  l”1^  |r2(u)f  in  ^(t.ujD.A),  shown  in  figure  3-2,  and 
x(t)  is  the  output  of  9  (t,u;D,A),  shown  in  figure  3-4.  Equation  (4-2), 
therefore,  provides  a  means  for  eliminating  the  operator  ^  |  •  J  from 

the  system  in  figure  3-4. 

Alternative  estimator  structures  can  be  obtained  by  straightforward 
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generalizations  of  the  material  In  reference  19,  pp.  15-23. 

Canonical  Realization  No.  2,  shown  In  figure  4-2,  Is  a  vector 
f 1 1ter-correlator  receiver.  The  matrix  Impulse  response  H'(t,u;D,7C)  is 
defined  by 

,  (H.(t,u;D,A)  ;  t  >  u 

H  (t,u;0,A)  =  <  '  (4-3) 

(0  ;  t  <  u  . 

Note  that  the  output  of  the  realizable  filter  H'(t,u;D,ft)  is  not  the 
causal  MMSE  estimate  of  s(t;d,‘3'),  given  d  =  D  and  Tt  = 


Canonical  Realization  No.  3,  shown  in  figures  4-3 ( a )  and  4-3(b),  are 
vector  versions  of  the  filter  squarer  receivers  shown  in  figures  2-5  and  2-6 
of  reference  19.  The  matrix  impulse  responses  Hfn(t,u;D,fr)  and 
Hf-c(t,u;0,^')  are  the  noncausal  and  causal  solutions,  respectively,  to 


Hn(t,u;0,A) 


T. 


Hf(Z,t;0,A)Hf(Z,u;0,A)dz, 


Ti  < 


<  Tf 


(4-4) 


As  with  the  scalar  case,  there  are  an  infinite  number  of  noncausal 
solutions  because,  with  H  (t.ujD.’/C)  given  by  (2-20), 


Hfn(t,u;D,A) 


X.(0,A) 

X^O.A)  +  N  q/2 


0.(t;D,A)a.  ( u ; 0 , A) ,  T.  <  t,  u  <  Tf 


(4-5) 


is  a  solution  to  (4-4)  for  any  assignment  of  plus  and  minus  signs.  The 
substitution  of  (4-4)  into  (2-19)  yields 
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5.0  THE  SPECIAL  CASE  OF  CONSTANT  DELAY,  STATIONARY  PROCESSES,  AND  LONG 

OBSERVATION  INTERVAL  (CDSPLOT) 

This  section  links  the  general  solution  to  the  problem  of  ML  time-delay 
estimation  presented  above  to  the  solution  of  Knapp  and  Carter1  in  which 
the  delay  is  constant,  d(t)  =  dQ,  the  signal  process  is  stationary,  and 
the  observation  interval  is  long.  We  call  this  the  CDSPLOT  case.  This 
exercise  has  the  benefit  of  providing  additional  insight  into  Knapp  and 
Carter's  solution,  as  well  as  a  more  explicit  description  of  the  bias  term. 

As  was  shown  in  (2-16)  the  log-likelihood  function,  lnA(D,A*),  consists 
of  the  sum  of  a  data-dependent  term,  lR(D,/C),  and  a  bias  term,  fB(D,A). 
The  forms  of  these  terms  under  the  CDSPLOT  approximation  are  derived  in  the 
next  two  subsections. 

5.1  DA TA-OE PENDENT  TERM  i  (D,A)  UNDER  CDSPLOT  APPROXIMATION 

The  data-dependent  term,  lR(D,A),  is  given  in  the  general  case  in 

(2-19),  where  the  entries  of  H  (t,v;D,/C)  are  given  in  (3-49),  (3-50), 

— n  — 

(3-53),  and  (3-54).  It  can  be  seen  that  the  entries  themselves  are 
specified  in  terms  of  9n(t,u)  and  k(u)  in  (3-44)  and  (3-31), 

respectively.  Under  the  CDSPLOT  approximation,  we  obtain  gn(t,u) 
(approximately)  by  replacing  T.  and  T  in  (3-44)  by 
respectively,  R  (t,u)  by  R  (t  -  u),  and  Q(u)  in  (3.41)  by 


- «  and  +oo, 
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I 

Q(u)  =  -§ 

where 


I 


A 


0 _ 

+  1 


-  00  <  u  <  00  , 


(5-la) 


(5-lb) 


(Equation  (5-1)  is  obtained  by  setting  D(t)  =  DQ  in  (3-41)  and  noting 
that,  with  [T..,Tf]  =  [-*>,  »  ],  the  middle  expression  for  Q(u)  in  (3-41) 
applies  for  all  -*><  u  <».)  With  (3-44)  so  modified,  we  try  a  solution  of 
the  form  9n(t,o;0,A)  =  gn(t  -  o;D0,y).  Thus,  under  the  CDSPLOT 
approximation,  (3-44)  becomes 


+  « 


Rs(t  -  u) 


=  /9" 


(t  -  <j;D0,A) 


Rs(o  -  u) 


do 


*  2  9n(t  -  u;D0.A)  (5_2a) 


for  -  *><  t ,u  < *>  . 


The  above  can  be  notationally  simplified  by  a  change  of  variables  and 
by  again  assuming  the  dependence  of  9n(*)  on  Oq  and  T  implicitly.  This 
leads  to 


Rs(t)  « 


*  /  9n<  T  - 


X) 


Rs00  dx  +  gn(r) 


for  -  oo<  T  <  oo  . 


( 5 -2b ) 


The  solution  to  ( 5 -2b )  can  be  obtained  easily  by  Fourier  transforms  and  is 


ir*  .V  **  /*» /** .V.VAVimS  AW^^ / »  /(« 


i  »■*  h*  lb*  ■■  IT  U"  V .!» . 
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9p(t) 


G(f)  e^2*ftdf  * 


(5-3) 


where 


G(f)  = 


S$(f) 

S$(f)  +  N^/2 


(5-4) 


and  where  S$(f)  is  the  power  spectral  density  of  s(t).  Since  s(t)  is 

real,  S  (f)  and  G(f)  are  real  and  even.  This  implies  that  g  (t)  is  real 
s  n 

and  even: 


y-1)  -  yt>. 


(5-5) 


It  can  also  be  seen  from  (5-2)  that  9n(t)  does  not  depend  on  Dq. 
Also,  with  D ( t )  =  Dq  and  with  T..  *  Tf  *  +•»,  k(u)  in  (3-31)  reduces 
to  the  constant 


A2  -  1 

k(u)  =  zz - 1  «  k- 

A^  ♦  1 


(5-6) 


The  substitution  of  the  above  with  [T.,Tf]  =  [-«,  ®  ],  into  (3-49), 

(3-50),  and  (3-53)  will  cause  (3-54)  to  yield  the  t ime-i nvariant  impulse 
responses , 


hn(t)  -  p— Vt} 

A  +  1 


h12(t)  =JT~9n(t  *Do} 

A  *  1 


(5-7) 


(5-8) 


>  ^vsv-v^v >y vVrVvivi'I.Xvivf-' 
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*21™  "  i*T7  9n(t  -  °0> 


h„(t)  =  7?  9n0:)  * 

“  tc  +  1  n 


which,  when  substituted  into  (2-19),  yield 


(5-9) 


(5-10) 


+  oo  +oe 


*R  (Do*A)  =  t 


-  v)  R^t)  R1  ( v)  dtd v 


4*«o  4*00 


+  r0  j  j  PT7  V1  -  *  ♦  Do>  Vl>  "i")  «d* 


+  oo  +oo 


* J  J  JT77  V*  -  V  ¥‘>  ¥*> dt'"' 


+  00  +00 


r0jj  ^77  V1  - v)  ¥*>  Vv> dtdv 


(5-11) 


In  (5-11)  we  have  written  the  integration  limits  as  + »  for 
convenience.  Since  r(t)  is  defined  as  zero  outside  [T.,T^],  the 
integration  is  actually  still  over  the  long  but  finite  interval 
[T.,T^].  Therefore,  the  integrals  exist. 


We  are  interested  in  finding  the  ML  estimate  of  d  ,  D  ,  which 

o  oJMl 

will  be  obtained  by  choosing  D  to  maximize  1_(D  ,/C)  +  l_(D  , A) . 

o  R  o  B  o 

We  will  find  in  the  next  subsection  that  lQ(D„,/C)  does  not  depend  on 

D  0 
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iD(0  ,/C)  In  (5-11).  Note  that  only  the  two  middle  terms  in  (5-11) 

K  0 

depend  on  0Q.  Using  the  fact  that  9n(t)  is  an  even  function,  we  note 
that  these  two  terms  can  be  combined  and  written  as 


+  oo 

ip  (0o,A)  =  jjj  jj  gn(t  -  a)  R1  (t)  R2(d  +  D0)  dtdv  .  (5-12) 

-  oo 

Equation  (5-12)  can  be  written  another  way  by  introducing  functions 
h.j(t)  and  h2(t),  satisfying  the  equation 

2  A 


The 

substitution  of  (5-13) 

into  (5-12)  yields 

+  oo  4-  oo 

C  C 

4*  OO 

r 

JR  <V> 

=  I  dz  1  h] ( z  -  t) 

R-|(t)  dt  h2(z  - 

a)  R 2(a  +  Dq)  do  ,  (5-14) 

=  OO  -OO 

-  oo 

and  we 

see  that  l.(D  ,A“) 

R  0 

can  be  obtained 

from  the  "generalized 

correlator"  shown  in  reference  1,  figure  1.  By  taking  the  Fourier  transform 
in  (5-13)  it  follows  that 

iT  'TA~~  G<f>  =  Mf)  Mf)  E  +<f>  •  (5-15) 

o  Jr  +  l  1  c 

where  +(f)  is  the  "frequency  weighting  function"  appearing  in  reference  1, 
equation  (6).  Combining  (5-4)  and  (5-15)  gives 


g„(t  -  .) 


■: 


h-,  ( z 


t)  h2(z 


-  o)  dz 


(5-13) 
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♦(«  ■  S2  Jr4 


ss(f> 


No  A2  +  1  S  If)  +  NV2 
S  0 


(5-16) 


The  function  (f)  can  be  written  in  terms  of  the  coherence  function  of 
r^t)  and  r2(t),  defined  as 


S  (f) 
r  r  v  7 

...  a  rr  2 _ 

Ti2(f>  ~  ( - 

^sr  (f>  Sr$(f) 


(5-17) 


where 


S  (f)  =  power  density  spectrum  of  r..(t), 
1  1 


Sr/f>  -  Vf>  *1  ■ 


(5-18) 


S  (f)  *  power  density  spectrum  of  r  (t), 
r2  ^ 


S_  (f)  *  A2  S  (f)  ♦  -§  . 

2  S  * 


(5-19) 


S  (f)  =  cross-power-density  spectrum  of  r.(t)  and  r„(t), 
12  1  l 


S_  (f)  =  A  S,(f)  e 
12  s 


+j2irf  0 


(5-20) 


The  detailed  algebraic  steps  are  shown  below,  which  starts  by 
substituting  (5-2)  into  (5-16): 
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•Kf)  = 


A  _2 
+  S2  No 


ss(f) 


Ss(f)  + — ^r 
S  2(1  +  A  ; 


AS$(f) 

9  N„ 

(1  +  h  )  -I  S  (f ) 


ASs(f) 


N  "  "  N  " 

|ss(f)  *  -fj  [a2  Ss(f)  +  -f  - 


sr  r  (f)  SGN(A) 

r1  r2  I  _ 

S_  (f)  S  (f)  -  S2  (f) 

rl  r2  rl  r2 


A2  S$(f) 


(5-21) 


By  multiplying  the  numerator  and  denominator  in  (5-21)  by 


Sr  r  (f)  /[Sr  (f)  Sr  (f)] 
rTr2  r]  r2 


we  obtain 


*(f)  = 


|t12(f)|  SGN(A) 


|Sr,r2<f>|  f  -  h2<f>|2l 


(5-22) 


For  A*  >  0,  this  is  the  frequency  weighting  function  associated  with  the  ML 

2  3 

or  HT  (for  Hannon/Thomson  ’  processor  in  reference  1,  table  1.  Knapp  and 
Carter1  do  not  include  the  factor  SGN(A)  because  their  receiver  has  a 
square  law  device  before  the  peak  detector. 
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5.2  BIAS  TERM  lD( 0  .£)  UNDER  CDSPLOT  APPROXIMATION 
o  0 


The  bias  term  for  the  general  case  Is  given  by  (3-56),  where  gc(t,u) 
Is  the  solution  to  (3-44)  for  Tf  =  t,  with  gjt.u)  =  0  for  t  <  u.  Under 
the  CDSPLOT  approximation,  we  set  0(t)  =  Dq,  R$(t,u)  =  R$(t  -  u),  and 
Tf  *  t  in  (3-44)  to  obtain 


[t  -  u)  =  j  gc(t. 


a)  R$(o  -  u)  do  +  Q(u)  gc(t,u) 


(5-23) 


T.-D 
i  o 


for  T.-D  <  u  <  t.  Under  the  CDSPLOT  approximation  the  function  Q(u) 
1  o  - 

in  (3-41)  becomes 


Q(u)  = 


-f  :  Ti  -  do  <  u  <  Ti 


;  T.  <  u  <  t  -  D0 


2  ;  t  -  0Q  <  u  <  t  , 


(5.24) 


where  Nq  was  defined  in  (5-lb).  Recall  that  gc(t,o)  is  the  impulse 

response  of  the  casual  LMMSE  estimator  of  s(t)  from  z(t)  =  s(t)  f  n(t), 

where  n(t)  has  covariance  function  Q(u)6(t  -  u)  (see  (3-36)  through 

(3-40)).  Looking  at  (5-24)  we  see  that  under  the  CDSPLOT  approximation  n(u) 

is  "piecewise"  stationary  in  the  three  intervals  [T.  -  Dq,  T.], 

[T.,  t  -  Dq],  and  [t  -  Dq,  t];  but  0(u)  changes  abruptly  at  the 

interval  boundaries.  If  we  let  T.  -»  -«,  then  n(u)  will  be  stationary  for 

-«<  u  <  t  -  0  and  g  (t.u)  will  be  time-invariant  in  this  range.  Thus, 

-  o  c 
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9c(t,u)  =  g^(t  -  u)  for  -  *  <  u  <  t  -  0Q  , 


(5-25) 


t 


where  the  function  '(v)  is  the  solution  to  the  Wiener-Hopf  equation: 


Rs(t) 


00 

-J* 


(v)  Rz( r  -  v)dv  ;  0  <  t , 


(5-26) 


R,(t)  =  Rc(r)  +  -f  6(t) 


(5-27) 


3 


g  j  (T)  =  o  ;  r  <  0  . 


(5-28) 


Since  the  statistics  of  n(u)  change  abruptly  at  u  =  t  -  D  and 
gc(t,u)  operates,  in  general,  over  all  past  data,  we  cannot  expect 
gc(t,u)  to  be  time-invariant  for  t  -  D  <  u  <  t.  Thus,  under  the 
COSPLOT  approximation,  g  (t,u)  is  a  time-varying  casual  impulse  response 
that  has  the  approximate  time-invariant  form  g^(t  -  u)  specified  in 
(5-26)  for  T.  <  u  <  t  -  D0>  but  not  otherwise.  Using  these  results  in 
(3-56),  with  f(t)  *  t  -  0  ,  we  have 


V°o*A>  ~  ~  2 


-  1 1 


Vo 


o)do  . 


rTi  .  v 

•'ft  r  » 


(0)do 


(5-29) 


Because  0  is  finite,  the  values  of  the  first  and  the  third  integrals  in 
0 

(5-29)  are  negligible  compared  with  that  of  the  second  integral  for 


m 
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sufficiently  large  t  -  .  Consequently, 


16  <0o.A)  =  -  i  ct  -  0o  -  T, ]  g'11  (0) 


*  -  *  [t  -  t,j  g<’>  (0)  . 


(5-30) 


We  can  obtain  a  more  explicit  form  for  iB(D  .ft)  by  referring  to 

D  0 

(3-57),  which,  under  the  COSPLOT  approximation  (for  =  -  *>  <  u  <  t  -  0Q) , 
becomes 


eoc(u)  =  E{( s( u) —  sc(u))2} 


-i«>> 

=  (5-3i) 


If  the  signal  spectrum  ( f )  Is  rational  with  finite  variance,  then 
(reference  17,  p.  501) 


+  9 

■>[ 


In  [1  +  S-  S_(f ) ]  df  . 


(5-32) 


Combining  (5-30),  (5-31),  and  (5-32),  we  have 


*8  <Do 


V 

.A)  =r  -  1  Ct  -  Ti]  I 


In  [1  +  ,  S,(f)]  df  , 


(5-33) 


and  we  see  that,  as  in  reference  1,  the  bias  does  not  depend  upon  Dc 
for  COSPLOT. 


1 T 
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6.0  SUMMARY 

This  report  has  generalized  previous  theory  concerning  ML  time-delay 
estimation  to  include  time-varying  delay,  finite  observation  interval,  and 
nonstationary  signal  process.  It  has  presented  several  receiver  structures 
that  can  be  used  to  obtain  the  ML  estimates  of  time-delay  and  attenuation  in 
one  of  two  received  signals  compared .  with  the  other.  Here  it  is  shown  that 
the  general  theory  reduces  to  that  of  Knapp  and  Carter1  for  constant 
delay,  stationary  signal  process,  and  long  observation  interval. 
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